1_Dynamics/2_Time_domain_simulation/Fed batch bioreactor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fed Batch Bioreactor\n",
    "This model represents the fed batch bioreactor in section 2.4.9 of Seborg et al."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Nr  | Symbol      | Name                                 |                Units|\n",
    "|-----|-------------|--------------------------------------|---------------------|\n",
    "| 1   | $\\mu_\\text{max}$ |   Maximum specific growth rate  |                1/hr |\n",
    "| 2   | $K_s$       |   Monod constant                     |                 g/L |\n",
    "| 3   | $Y_{X/S}$   |   Cell yield coefficient             |                   1 |\n",
    "| 4   | $Y_{P/S}$   |   Product yield coefficient          |                   1 |\n",
    "| 5   | $X$         |   Cell mass concentration            |                 g/L |\n",
    "| 6   | $S$         |   Substrate mass concentration       |                 g/L |\n",
    "| 7   | $V$         |   Volume                             |                   L |\n",
    "| 8   | $P$         |   Substrate mass concentration       |                 g/L |\n",
    "| 9   | $F$         |   Feed flow rate                     |                L/hr |\n",
    "| 10  | $S_f$       |   Substrate mass concentration       |                 g/L |\n",
    "| 11  | $\\mu$       |   Specific growth rate               |                1/hr |\n",
    "| 12  | $r_g$       |   Cell growth rate                   |        g /( L hr)   |\n",
    "| 13  | $r_p$       |   Product growth rate                |          g /( L hr) |  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Model equations:\n",
    "    \n",
    "| Nr     | Equation                              | Inputs     | Outputs         | Parameters              |\n",
    "|--------|---------------------------------------|------------|-----------------|-------------------------|\n",
    "| 1      | $$r_g = \\mu X$$                       |      -     |        $r$, $\\mu$, $X$ |   -  |                    \n",
    "| 2      | $$\\mu = \\mu_{max} \\frac{S}{K_S + S}$$ |      -      | $S$             | $\\mu_\\text{max}$, $K_s$ |\n",
    "| 3      | $$r_p = Y_{P/X}r_G$$                  |      -      | $r_p$           | $Y_{P/X}$               |\n",
    "| 4      | $$\\frac{d(XV)}{dt} = Vr_g$$           |      -      | $V$             |   -                      |\n",
    "| 5      | $$\\frac{d(PV)}{dt} = Vr_p$$           |      -      | $P$             |    -                     |\n",
    "| 6      | $$\\frac{d(SV)}{dt} =  FS_f - Vr$$     | $F$, $S_f$ |  -               | $Y_{X/S}$               |\n",
    "| 7      | $$\\frac{dV}{dt} = F$$                 |      -      |  -               |    -                     |\n",
    "| Number | 7                                     | 2          | 7               | 4                       |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_max = 0.2  # Maximum growth rate\n",
    "K_s = 1.0  # Monod constant\n",
    "Y_xs = 0.5 # Cell yield coefficient\n",
    "Y_px = 0.2  # Product yield coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "States"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = 0.05  # Concentration of the cells\n",
    "S = 10   # Concentration of the substrate\n",
    "P = 0   # Concentration of the product\n",
    "V = 1   # Reactor volume\n",
    "\n",
    "x0 = [X, S, P, V]  # State vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = 0.05  # Feed rate\n",
    "S_f = 10  # Concentration of substrate in feed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dxdt(t, x):\n",
    "    [X, S, P, V] = x\n",
    "    \n",
    "    mu = mu_max * S / (K_s + S)\n",
    "    rg = mu * X\n",
    "    rp = Y_px * rg\n",
    "    dVdt = F\n",
    "    dXdt = 1/V*(V * rg - dVdt*X)\n",
    "    dPdt = 1/V*(V * rp - dVdt*P)\n",
    "    dSdt = 1/V*(F * S_f - 1 / Y_xs * V * rg - dVdt*S)\n",
    "\n",
    "    return [dXdt, dSdt, dPdt, dVdt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tspan = [0, 30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsmooth = numpy.linspace(0, 30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fs = [0.05, 0.02]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "for F in Fs:\n",
    "    out = scipy.integrate.solve_ivp(dxdt, tspan, x0, t_eval=tsmooth)\n",
    "    results.append(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['X', 'S', 'P', 'V']\n",
    "units = {'X': 'g/L', 'S': 'g/L', 'P': 'g/L', 'V': 'L'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hVVdbA4d9K7wkklFBCqNJ7FcWuiL2j2HVwRsc2ttGZsY0zOvM5Y2/YK/ZRFBBBaSq9dwgdEkJIQnrP+v44F0VMuQnJLcl6n+c+uWWfc9cJ2axz9tlFVBVjjDHG1wR4OwBjjDGmKpagjDHG+CRLUMYYY3ySJShjjDE+yRKUMcYYn2QJyhhjjE+yBGWMjxKRN0Rkv4isreZzEZFnRSRFRFaLyGBPx2hMYwrydgDuSEhI0OTkZG+HYQwAy5YtO6CqrTzwVW8BzwPvVPP5mUB312ME8JLrZ42sPhlfU12d8osElZyczNKlS70dhjEAiMhOT3yPqs4TkeQaipwHvKPOaPuFIhInIomqmlbTfq0+GV9TXZ2yJj5j/Fd7YPdhr/e43vsNEZkoIktFZGlGRoZHgjPmaFmCMuZwqrD2M0hd4e1I3CFVvFfl3GWqOklVh6rq0FatPNE6aYxL+jpY9VG9NrUEZcwhWdvh/Yvh0+th8avejsYde4COh73uAKR6KRZjfq20AL79G7x8PHz3KJSX1HkXfnEPyphGpQoLnofvH4OAIBj7BAz7nbejcscU4I8i8iFO54ic2u4/GeMR2+fDFzdDzi4YdCWc9ncICq3zbixBGTP3XzDncThmHIx7EmKrvI3jcSIyGTgRSBCRPcBDQDCAqr4MTAPGASlAIXCddyI15jDb5zstEbEd4NppkDy63ruyBGWat5+ed5LTwAlw7vMQ4Dut3qp6eS2fK3CLh8IxpnZ7lsLk8RDXCa6bBpEJR7U736mNxnjasrfh279A7/PgnGd9KjkZ43f2rYX3LnKS0tVfHnVyAktQprnavRi+vgO6nQYXvgaB1phgTL0V58LkyyE4Aq6eAjGJDbJbq5Wm+SnJh88nQkwHuPgNCArxdkTG+Lfp90HuHrh+BrTo1GC7tQRlmp9v/wLZO5w28rAYb0djjH9b/yWs+gDG3Asdhzforq2JzzQvm6bDsrdg9O3Q6VhvR2OMf8tNg69uh3aD4IR7G3z3lqBM81GYBVNuhTb94KQHvB2NMf5N1UlOZcVw4asQGNzgX2FNfKb5+P4xJ0ld9b96DRo0xhxm49ewZQac/g9I6N4oX2FXUKZ5SF0BS9+A4ROhbT9vR2OMfystgG/uh9Z9YMTvG+1r7ArKNH2VlTD1LohsBSfd7+1ojPF/856EnN1w3fRGHaJhCco0fSvehb3L4IJJEBbr7WiM8W8HtsBPz8GAyxu9o5HXmvhEJFBEVojI196KwTQDhVkw62FIOhb6X+rtaIzxb6ow/V5nQO5pjzb613nzHtTtwAYvfr9pDub+C4oPwrj/A6lq+SRjjNs2fwNbv3d6wUa1dmuTsopKnGkj684rCUpEOgBnAa954/tNM5GxGZa8BoOvgbZ93d7s8ekbmLk+vREDM8YPlZfCjL9AQg8YdoPbmz0zawtXvLqIkvKKOn+lt66gngbuBSqrK2BLVJuj9u1fnaaIk/7i9iYLtmbyytxtrN5zsBEDM8YPLXkVsrbCGf90e8zTtox8Js3bRtvYMEKDAuv8lR5PUCJyNrBfVZfVVM6WqDZHJWWWM0ZjzN0Q5d7fT1lFJQ9+uZYOLcK5+cRujRygMX6kIBPm/Au6nQrdT3NrE1XloSnrCA0O4P5xPev1td64ghoNnCsiO4APgZNF5D0vxGGaqopypymiRXKdxmi8+eN2tuzP56Fz+hAeUvezPWOarNn/gNJ85+rJTdPW7GP+lgPcffoxtI4Oq9fXejxBqer9qtpBVZOB8cD3qnqlp+MwTdiyNyFjY52WmU7LKeLpWVs4pWdrTuvdppEDNMaPpK936tSwG6DVMW5tkl9Szt+/Xk+fdjFcObL+s5vbOCjTtBRlw+x/QvLx0Osctzd7bOoGKiqVh8/t04jBGeNnVGHG/RAaAye6P8j9mVmb2ZdbzItXDiYwoP69Z7061ZGqzlHVs70Zg2li5j3pJKkz/ul2t/L5WzKYujqNW07qRseWEY0coDF+ZPM3sG2O0608oqVbm2zcl8sbP+5g/LCODE5qcVRfb3PxmabjQAosehkGXwWJ/d3apLisgr99sZbOCZFMHNOlkQM0xo8c3q186PVubVJZqfz1f2uJDQ/mvrH16xhxOGviM03HzL9BUDic/De3N5k0bxs7Mgt55/rhhAVbxwhjfnaoW/mET93uVv7Z8j0s3ZnNvy/uT4vIo1+p2q6gTNOw9XvYNA3G3OX2CPedmQU8PzuFs/onMqaHDWUw5mf5GXXuVp5dUMrj0zcytFMLLh7coUHCsARl/F9FGUy/z9Wt/A9ubaKqPDxlHcEBwt/O6t248Rnjb757BMoK4IzH3d7k3zM2klNUxt/P70vAUXSMOJwlKOP/Fk+CA5th7BMQ7N54i2lr9jF7UwZ3ntaDtrH1G6NhTJO0dxmseA9G/gFa9XBrk8Xbs5i8eDc3HNeZXokxDRaKJSjj3/L3w5wnoNtp0GOsW5vkFJXx8Ffr6Ns+hmuPTW7c+IzxJ5WVMO0ep5l8zL1ubVJSXsED/1tD+7hw7ji1YVfWPapOEiLSAmgHFAE7VLXaufWMaRSzHoGyIufqyc1u5f/+ZiOZ+SW8cc0wggLtHM2Yn636wLV22isQ5t6V0MtztpGyP583rxtGREjD9rur895EJBa4BbgcCAEygDCgjYgsBF5U1dkNGqUxVdm9BFa+B6NvhwT35s5btjOb9xft4vrRnenXwTOLF4pIADCAX07m1qmqTZdufEtRtrN2Wofh0P8ytzbZmpHPC7NTOGdAO046xr3OSXVRn3T3KfAOcLyq/mrKZxEZAlwlIl1U9fWGCNCYKlWUw9d3QnQ7GHOPW5uUllfywOdraBcbxl2nu9e2fjREpCtwH3AqsIVfTuZ6iEgh8ArwtrU8GJ8w6xEozIQrP3OrNaKyUrn/8zWEBQfw4NmN09GozglKVWvqc5imqnccRTzGuGfRy5C+Bi59F0Kj3drkhdkpbErP47WrhxIZ6pEhgI8BLwE36RErtolIa+AK4CrgbU8EY0y1di925tsbeQskDnBrk/cW7WTx9iz+fVF/WkW7N+dlXTV0LV0IJDXwPo35tYO7nfn2up/h9nx761NzeWF2CucPbMepHpoMVlUvr+Hjzqr6tEcCMaYmFWXw1R0Q096Z0sgNuzILeWL6Rsb0aMUlQxtmzFNVGvoOsa2pbRrfN38GrXR7Gfeyikru/mQVcREhPHSOz0wG+4m3AzAGgIUvwv51cOa/ITSq1uKVlcp9n60mQIQnLuyHuNk5qT4aOkHVb+F5Y9y14WvY+DWceB+0cG8a/5fnbGV9Wi6Pnd+nQaZfaSB2Mme8L2u7M0zjmHHQy715uz9YvIsF2zL5y1m9aBcX3qjh1acX33NUnYgEiDvqiIypTmEWTP0TtOkHo/7o1ibrU3N59vstnN0/kbF9Exs5wDqxkznjXZWVMOVWkECnNcIN2w8U8I+pGziuWwLjh3Vs5ADrdw9qaT0/M+bozHjA6WU04RO3Jq8sLqvgjo9WEBcRwqPn9fVAgL8mIl9R/clcvIfDMebXlr0BO+bDOc9CbO33kcoqKrnjo5WEBAXwf5f0b9SmvUPqk6DaA9NVdUVDB2NMtTbPgFWTnS7lbvYy+tc3G9mcns/b1w+npXea9p6s52fGNK7snTDzIehyEgy+2q1NnvtuC6t2H+SFKwaTGNu4TXuH1CdBbQNuF5EBwCpgOvCtqmY3aGTGHFKc4/QyatXL7TFP8zZn8OaPO7j22GRO8N5M5RNw6scsVc3zVhDG/IoqfHWb8/zcZ93qaLR0RxbPz07hosEdOKu/55rK69xJQlU/VNVrVXUQ8AzQBfhcROaJyIMiMrzBozTN27R7IH8fnP8CBNU+3iKroJS7P1lFt9ZR/PnMo1807Si8gTODxDQR+U5E7nOd2LlNRMaKyCYRSRGRP1fx+bUikiEiK12PGxsqeNNELZ7krJJ72qMQV/uooJzCMu74aCXtW4Tz8Lmenfn/qMZBuZr5VgCPi0gMcBpwI7C4AWIzBtZ8Cqs/ghPvh/ZDai1eWanc8dFKDhaV8eZ1w7y6CKGqLsQZG/iwiMQDpwN3iUh/YDnwjap+XN32IhIIvIBTr/YAS0RkiqquP6LoR6rqXq8R07ylr4Nv/wbdT3drlVxV5e5PV7Evp5iPfz+K6DD3Fi5sKPVOUCJyYRVvHwT+Wv9wjDnMwV3w9Z+cucGOv9utTV6au5V5mzN47Py+9Gnnmbn23KGqmcBk1+PQtGC1Tb8+HEhR1W2ubT4EzgOOTFDG1K6sGD670ZkE9rwX3Wrae23+dmauT+dvZ/dmcFILDwT5a0dzBXUDMAo4NDHsiThniz1E5FFVffcoYzPNWWUFfH4TaAVcOAkCa/9TXbA1k/98u4lzB7RjwgjfmdBERP5Uxds5wNRaNm0P7D7s9R5gRBXlLhKRMcBm4E5V3X1kARGZCEwESErynd+N8aBZD8P+9c4S7lG135dduiOLJ77ZyNg+bbl+dHKjh1eVoxmoWwn0UtWLVPUioDdQglOB7muI4EwzNu9J2PWTMz6jZedai+/PK+a2D1eQHB/JPxt5dHs9DAV+j5Nw2uMkihOBV0WkpkV3qjqII7utfwUkq2p/YBbVzOunqpNUdaiqDm3Vypa3b3Y2TYdFL8Hwm9xawn1/XjF//GAF7ePC+beHupRX5WgSVPIRSwbsB3qoahZQdnRhmWZt6/cw53HodykMqGk6O0dJeQW/f3cZ+cXlvHjlYKI8MxFsXcQDg1X1LlW9CydhtQLGANfWsN0e4PDRkB2A1MMLqGqmqpa4Xr4K1H6jzjQvWdud1ojEAU7HiFocqk85RWW8dOVgYjx83+lwR1OT54vI1/wyp9hFwDwRicS5F2VM3eXscdrJW/WEc56utZ1cVfnr/9ayfNdBXpwwmJ5tG2656QaUBJQe9roM6KSqRSJSUs02AEuA7iLSGdgLjMeZAf1nIpKoqmmul+cCGxoubOP3yorg46uda/FL34HgsBqLH1mfvH0f92gS1C3AhcBxOIf/DvCZa1mBkxogNtPclJfCJ9dCeQlc9i6ERNa6yZs/7uCTZXu47eRujOvnU1MZHe4DYKGIfOl6fQ4w2XUyV22HB1UtF5E/AjOAQOANVV0nIo8CS1V1CnCbiJwLlANZ1HxFZpqbaffAvtVw+UfQIrnW4j/Xp1O6+0R9qs9cfKIuwGeuR5VlGiJA04x882fYswQufhMSutdafM6m/fxj2gZO792GO05t/AUI60tV/y4i0/jlZO73qnpoWrAJtWw7DZh2xHsPHvb8fuD+ho3YNAlLXocV78Lxd8ExtXUYhZnr03ls6nrO6NOGO06pvf55Qn2uoGaLyGfAl6q669CbIhKCUwGvwenZ91aDRGiah8WvwtLX4dhboW9VIxh+bc2eHG5+fznHtInmv5cNJCDApzpFACAiUaqaD6Cqy4BlNZUxpsFsne1cPXU/A076S63Fl+/K5tbJy+nXPpanfKg+1aeTxFigAqeJIlVE1ovIdpwlrS8HnlLVt6rbWEQ6ishsEdkgIutE5PZ6RW6ajpTvYPp90ONMOPWRWovvzirkureW0CIihLeuG+aLnSIO+VJE/iMiY1zNeQCISBcRuUFEZlD7WChj6ubAFvjkGmh1DFz0GgTUPFh9+4ECbnx7KW1iwnj92mFEhPhOfarPku/FwIvAiyISDCQARarqbseIcuAuVV0uItHAMhGZWcXoeNMcZGyGT65zOkVc9GqtlSmroJRr3lxMWUUlH04cQeuYmm/6epOqniIi44CbgNEi0hKng8QmnDFQ16jqPm/GaJqYwiz44DIICIbLP3QG5dYgPbeYa95wJv55+7rhJEQ1ztLt9XW0Ux2VAWm1Fvz1NmmHtlHVPBHZgDM2xBJUc5ObBu9d5CydcflkCI2usXhOURlXv7GIPdlFvHfDCLq1rrm8L6jqHpIxjaK0AD641OkJe/WXtS7oeSC/hCteXUhmfgnv/24kyQm1d0rytIZeUbdORCQZGAQsquKziSKyVESWZmRkeDo009iKsuG9C6EoC678tNbKlF9SzrVvLmbTvjxeuXIIwzu39FCgxviBijKnO/neZXDx69BpVI3FDxaWcuVri9h7sIg3rh3GwI6+udas1xKUiETh9AC8Q1Vzj/zcRr43YaWF8MF4yEyB8e9Du0E1Fi8qreDGt5ewek8Oz10+mJN6tvZQoMb4gcpK+OJmSJkFZz8Nvc6psXhOYRlXv7GYbQcKePXqoYzo4rtrZ9Y5QYlItev8isjxbu4jGCc5va+qn9c1BuPHyorh46tg9yJnjr0uJ9ZY/NCV06LtWfz30gGM7dvWI2Ea4xcqK2HqnbDmYzjlQRhyTY3FD+SXMP7VhWxMy+OlCYM5vrtvn/zX5wpqrojcKyI/378SkTYi8h7w39o2FmdSp9eBDapaa3nThJQVw0dXOmd65zwDfS6osXhOURlXvb6IpTuzefqygZw3sL2HAm0YIhImIneIyPMictPhdcaYo3YoOS17yxnrdFxVcxL/Yl9OMZe9soDtB/J57ZqhnNKrjWfiPAr1SVBDgK7AChE52dVNfDGwgKpnWj7SaOAq4OTDFlkbV484jD85dOWUMhPOebbWM71M1w3cdXtzeXHCYL9LTi5v48y7twY4E/iPd8MxTcaRyenkv9U4LdiOAwVc+soC0nNLeOf6EYzx3irTdVKfbubZwE2uxDQLZ/LKkaq6x83tf6DqWZpNU1WS7ySnrd87V061JKftBwq49s3F7Msp5tVrhnpzyfaj1VtV+wGIyOvYQp6mIVSUwZd/hNUfupWclu3M5nfvOBOXvHfjCJ/tEFGV+kx1FAf8C+dqaSwwDpguIrer6vcNHJ/xdwWZ8MElkLoCznsBBl1ZY/Hlu7K58W2nMk2eONIri6Q1oJ9n9XfNq+fNWExTUFroDMLd8q0zQ8SYe2pMTt+s3cftH64gMTaMt64b7pNdyWtSnzbx5TgDdW9R1XLgWxEZiDNwd6eq1r4+gmkecvbAuxdA9k647D3oeVaNxaeuTuNPH6+kbWwYb/thZarCABE51ENVgHDXawFUVX1y6nXjowqzYPJ4Z77Ks5+qccl2VeXFOVt58ttNDOwYx2tXDyXexwbhuqM+CWrMkc15qroSOFZEftcwYRm/t3cZTL4Cygrhqs8h+bhqi1ZUKv/5dhMvztnK4KQ4XvXTynQkVa15Wgxj3JWxGSZf5pz0XfIW9D6v2qIFJeXc/ckqpq/dx7kD2vHvi/sTFuyff4r1uQdV7b0mVX316MIxTcLaz5xxGZGtneTUpk+1RXMKy7j9oxXM2ZTB5cM78vC5fQgN8s/KZEyjSJkFn1wPQSFwzdeQVH1ftK0Z+dz83nK27M/jL+N6cePxnX1tdek6sW6vpuFUVjgr4c77P+g40mnWi6q+g8PyXdnc+sEK9ucV848L+jJhRM2zSRjTrKjCT8/BrIegdW9nOrC4pGqKKp8u28ODX64jLDiAt68f7vNjnNxhCco0jPz9zkq42+fCwCvh7P9CUNXNdJWVyivztvHkt5tIjA3j45tGMci/O0MY07CKsuF/f4DN06HXuXD+SxAaVWXR3OIyHvxiLV+sTGVkl5Y8fdkg2sb67iTKdWEJyhy9HT/ApzdA8UE493kYfFW1RXdnFXLvp6tZsC2Ts/ol8s8L+xEbHuzBYI3xcXuWwqfXOZMpj/0XjLip2p568zZncN9nq0nPLeZPp/XglpO6Eegjazk1BEtQpv7KS+D7x5xmiJZd4MrPoG3fKouqKu8v2sXj0zYgIvzron5cOrSjX7ePG9Ogykth3r9h/n8gpgNc/w10GFpl0dziMh6ftpHJi3fRtVUkn/3h2CbZCmEJytRP2mr4302wfz0MuQ5Of6zaJoiU/Xn89Yu1LNyWxXHdEvjXxf1pHxfu4YCN8WHp6+GL30PaKhg4AcY+DmGxvymmqkxZlcpjUzeQmV/CTWO6cOdpPfy2l15tLEGZuiktgLn/gp+eh8gEuOIT6HF6lUWLSit47vstvDp/GxEhQTx+YT/GD7OrJmN+VlroXDX99JyTkC57r9rZyDen5/HIV+v4MSWT/h1ief2aofTv4D+zQtSHJSjjvs0zYNrdcHAXDLoKTnsUIn67LlNlpfL5ir08OWMT+3KLuXhIB+4/s2eTGNtkTINQdWaDmHYPHNzpdCw67VGI/O3SF/vzinlq5hY+WrKLyNAg/n5eH64Y0alJ3WuqjiUoU7v9G2DGX2Drd5DQA66dBsmjf1NMVfkh5QCPT9vI+rRcBnSI5bkrBjEs2RYXNOZn6ethxgOwbTbEd3fGNnX+7UpFOYVlvP7DNl7/YTsl5ZVcPSqZ207pTsvIEC8E7R2WoEz1ctOcMU3L3oKQKDj9HzB8ojNg8DCqyoKtmTw1azNLdmTTPi6cZ8YP5Jz+7QhoBmd5xrjl4G6nPq14F0JjYOwTMPSG39SnnKIy3vxxO6//sJ284nLG9WvLPWf0pLP/T/1VZ5agzG/lZ8CPT8OS16Cy3Jnz68T7f9P8UFmpzN60n5fnbmXJjmzaxITy6Hl9uGxYR5sNwphDclNh/n9h+dvO6+ET4YT7ftM8nnqwiDd+2M7kxbsoKK3g9N5tuOPUHvRu13ynbLQEZX6RvRMWvOCc4ZUXQ//xcMK90LLzr4oVlVbw1apUXp2/jS3782kXG8bD5/Rm/PCkJtubyJg6O7AFfnwGVn0IqDOT//F3Q9wvi5KrKst3ZfPOgp1MXZ2GAuf0T2TimK7NOjEdYgmquVN1JnZd9DKs/dwZENj/MjjuTkjo/quiOw4U8N7CnXyybA85RWX0bBvN05cN5Kz+iQQH1mftS2OaGFXYNgcWT4JN053ZVIZcA8feCi2Sfy6WU1TGV6tS+WDRLtan5RIdFsTVo5K5/rhkOrSI8Fr4vsYSVHNVWgDr/geLX4W0lRASDSP/ACNvhthfVq8tKCln6po0Pl26h8U7sggKEM7o05arRnViROeW1mXcGHDWPVvzMSx5HTK3QEQCjLkbht/083yUZRWV/JhygM+X72XGun2UlFfSs200/7ygH+cPakdEiP13fCT7jTQnlZWwexGsfA/WfQGl+dCqJ4x7EgaMh9BoAErKK5i7KYOvVqcxa306RWUVdE6I5J4zjuHiIR1oE9M05vky5qiUlzqrRK9837laqiyD9kPhgknQ53wICqW8opIlWzOZuiaVaWv2kVVQSkxYEJcN68glQzrSt32MneTVwBJUU6fqjE5f97nThJezG4Ijoc8FMGgCJI0CEfKKy5i7OpVv16Uze+N+8krKaRkZwoWD23PBoPYM6dTCKpIxFWWwY77T+rB+ijP/ZESCM1/ewAnQpjcFJeX8uOkAM9enM2tDOtmFZYQFB3BqrzacO6AdY3q0snu1brIE1RRVlMGuBbBxGmycCjm7ICAIup4MJ/8Vep6NhkSSsj+fOfO3M3vTfpbsyKKsQomPDGFcv0TG9U/k2K7xdm/JmKJs2DrbuUraMgOKc5xhFz3Pgr4XUdn5RDZkFPPTpkzmfrWIxduzKK2oJDosiFN6tuaMPm0Z06MVkaH2321d2W+sKVCF7O2wba6zuNm2uVCaB4GhTlI64R70mLPYVRzGou1Z/PS/Lfy0NZP9eSUAHNMmmutHd+aUXm0Y0qlFsxihbky1Kspg73Jn6ZiU75wl1rUCwltCz3OoOOZMNkUOY/HuQhYvyWLhR/PIKigFoHvrKK4bncwJx7RiaKeWhATZCd7RsATlj1SdLqy7FjiP7fMh17XQcWxH6HcRpcknszZsMMv3lbF8QzZLpq8gw5WQEqJCGNU1gdFd4xnToxXtbOJW05yVFkLqcld9Wgg7F0BZAQCaOJCC4bexPnI48ws7sXxPHqtW5JBfshSA9nHhnHhMK0Z3TWB0t4Qmsw6Tr7AE5etUIT/duY+0d5lzZrd3GRRlOZ9HxFPWcTR7ev6OlYH9+CknnjVbc0lZkE955WrAqUTHdo1nWHJLhnduSffWUXY/yTRP5aWQscGZjX/vMuexf70zIB0oadGD1A7nsjKoP98VdWfhPuHA9hJACQzYSa/EaM4f1I6hnVoyrHNLm5W/kVmC8iUl+XBgszP3XcYGSF8H+9ZAQQYAKgEUxXYjLf4E1gX24ofSbszPjCVtVYlrB8XER2bQt30sp/RqzYAOcQxMiqN1tJ3VmWamssKZhDVjs1OX9m+A9PVoxkaksgyA0qBoUiN7sSHmUn4o6co3OUlkpkVCGgQHCl1bhTKmRwz92sfSr30sfdrFEh5inRs8yRKUp5XkObOBZ++ArG2QtR0yU9DMFCR378/FKgJCyAhLZmvgEFaHJ/FTQTuWlSZRWOQkm/DgQLq2jmRk12iOaes8erWNoU1MqF0dmeahvNRp2s7e6SSjrG1o5jYqDmwhIHs7ARUlPxfNDGzNNunAyopxrCzrxDrtxM7iNgQWBpIUH0H3dlFcPjCabq2j6JkYTZeEKLt/5AMsQTUUVae3T/5+yEuDvH2Ql0ZlbiplWbvRnD0E5u4huPTgrzbLlyj2SCKbKrqwqXw0W7Udm7Uju7Q1QaXBJLWMIKlNBF1bRnBaq0g6J0SSHB9J+7hwm4i1GRCRscAzQCDwmqo+ccTnocA7wBAgE7hMVXd4Os4GpeqM0ctLh/x0NC+N4qw9lGbtpSJnLwG5ewkpSCW8JANBf96slCB2aRt2VLZmq55KirZna2U79gYnERfbivYtwklqGcHQ+Aguio+gc0IUHVuEE2Q9VX2WVxJUbZXO68pLoSTX6U5afJDS/CxK87IozcukvOAAFQWZUJBFQNEBgoqzCCnJJqIsk0At/82u8jSCVI0nTeNJ1SHs1tbs0VakB7QmPzKJiNgE2sSEkRgbTmJsGOfEhdE+Lpz2LcJJiAy1JNSMiUgg8AJwGrAHWCIiU1R1/WHFbo5NTxMAACAASURBVACyVbWbiIwH/gVc5vloq1FRDiW5lBRkU5ybRXFeFqX5ByjLz6KyIAstzEQKswgsziKkJIvwsmwiy7MJ1V+ufgQIB9AQ0rQlqRpPqvYijePIDG5LQUQHyqM7ENSiA21iI2kbG0aX2HBGx4bRLi6cFhHB1qrgpzyeoNysdG7Jy8kiNzOd8tIiKkuLKC8toqK0iIrSYirLCtHSIipLC6GsCMoKobQQKS8koKyAgLIiAisKCS4vILiikNCKQkIrC4ioLCCEsl99T4jr8fP3ajgHNYpMYsjUGLK0NQcDW1AYHE9RaDyl4a2pjGqLRLclOjqWlpEhxEeF0iUyhBHRobSKDiU23CqNqdVwIEVVtwGIyIfAecDhdeU84GHX80+B50VEVFWpo7LSEjJSd1BWUkRZcQHlJYWUlxRSUVLo1KuSArS00Onh5voZUFZIYHkBgeWFBJcXElJRQGhlIeGVhURoAeE4iSbU9ThyEfMiDSGTGLI1ihyJIS/wGApCWlASEk9peCsqItsg0W0JjmtPZExLWkaH0ioqlK5RIcRHhlozXBPnjSsodyqdW9ZNeZqRW59xu3yZBlJEKAWEUUwoxRJOjoRTEhhHSUA7yoIjKQ2KoiI4kvKQGDQ0Fg2LQcJbEBAeR1BkS0Ji4omKiCQ6LIi4sGCSwoOJDguyAa2mMbQHdh/2eg8woroyqlouIjlAPHCgrl+2N2U1yR+f6nb5Qg2liFCKJJzigHBKJILcwBhKQhIpD4qgPDiaypBoKkNj0LBYAsLiCIxsQVBkC0KjEwiLiScqKprY8GDaWB0yVfBGgnKn0rmlzZBzWBzdhoDgMCQ4jICgUAJCwgl0PYLCIgkOiyQoNILQ8GhCw0IJCwokOlDs6sX4g6r+SI+8MnKnDCIyEZgIkJSUVOWXtWzXhcUDHiMgJIyAYKcOBYdGEBQWRXBoBCHhUYSERxISHkVEZBQRQUHYvNumMXkjQTVYhercexidew9r0OCM8SF7gI6Hve4ApFZTZo+IBOG0omUduSNVnQRMAhg6dGiVzX8xcfEMv+DWBgjbmIbhjWtqdyodqjpJVYeq6tBWrVp5LDhjfMgSoLuIdBaREGA8MOWIMlOAa1zPLwa+r8/9J2N8kTcSlDuVzphmT1XLgT8CM4ANwMequk5EHhWRc13FXgfiRSQF+BPwZ+9Ea0zDE2+cbInIOOBpnG7mb6jqP2opnwHsrObjBOpxQ9gH2XH4lpqOo5Oq+u1lfS31CZrGv2FTOAZoPsdRZZ3ySoJqSCKyVFWHejuOo2XH4VuaynHUR1M49qZwDGDHYf06jTHG+CRLUMYYY3xSU0hQk7wdQAOx4/AtTeU46qMpHHtTOAZo5sfh9/egjDHGNE1N4QrKGGNME2QJyhhjjE/y6wQlImNFZJOIpIiI3wxQFJE3RGS/iKw97L2WIjJTRLa4frbwZoy1EZGOIjJbRDaIyDoRud31vr8dR5iILBaRVa7jeMT1fmcRWeQ6jo9cg8qbNKtP3mV16rf8NkEdtmzHmUBv4HIR6e3dqNz2FjD2iPf+DHynqt2B7/D9GQHKgbtUtRcwErjF9fv3t+MoAU5W1QHAQGCsiIzEWVfpKddxZOOsu9RkWX3yCVanjuC3CYrDlu1Q1VLg0LIdPk9V5/HbCT3PA952PX8bON+jQdWRqqap6nLX8zycqXja43/Hoaqa73oZ7HoocDLO+krgB8fRAKw+eZnVqd/y5wRV1bId7b0US0Noo6pp4PyhAq29HI/bRCQZGAQswg+PQ0QCRWQlsB+YCWwFDrrmwgP//9tyh9UnH2J1yuHPCcqtZTtM4xKRKOAz4A5VzfV2PPWhqhWqOhBnZv3hQK+qink2Ko+z+uQjrE79wp8TlFvLdviRdBFJBHD93O/leGolIsE4Fel9Vf3c9bbfHcchqnoQmIPT/h/nWl8J/P9vyx1Wn3yA1alf8+cE1dSW7Th8XZ9rgC+9GEutxFmS+HVgg6r+97CP/O04WolInOt5OHAqTtv/bJz1lcAPjqMBWH3yMqtTVVBVv30A44DNOO2bf/F2PHWIezKQBpThnLneAMTj9NDZ4vrZ0ttx1nIMx+Fcoq8GVroe4/zwOPoDK1zHsRZ40PV+F2AxkAJ8AoR6O1YP/C6sPnn3OKxOHfGwqY6MMcb4JH9u4jPGGNOEWYIyxhjjkyxBGWOM8UmWoIwxxvgkS1DGGGN8kiUoHyMi8SKy0vXYJyJ7D3v9UyN95yARec31/GERudvN7Wb5+szKpnmz+uTfgmovYjxJVTNxZgBGRB4G8lX1yUb+2geAx9wt7BpQKMC7wM3APxopLmOOitUn/2ZXUH5ERPJdP08Ukbki8rGIbBaRJ0RkgmsNljUi0tVVrpWIfCYiS1yP0VXsMxror6qrDnu7t4jMEZFtInKbq1yya52aF4HlONPiTAEub+zjNqYxWH3yfXYF5b8G4EzAmAVsA15T1eHiLHJ2K3AH8AzO+is/iEgSMIPfTto4FGe09+F6AicB0cAmEXnJ9f4xwHWqevOhgiISKiLxrjNVY/yV1ScfZAnKfy1R1xT8IrIV+Nb1/hqcygDOHFi9nRYEAGJEJFqdtWYOSQQyjtj3VFUtAUpEZD/QxvX+TlVdeETZ/UA7wCqU8WdWn3yQJSj/VXLY88rDXlfyy79rADBKVYtq2E8REFbDvisO219BFduHufZhjD+z+uSD7B5U0/Yt8MdDL0RkYBVlNgDd6rNz183dtsCO+mxvjJ+x+uRhlqCattuAoSKyWkTWA78/soCqbgRiXTd362oIsFB/WSXTmKbM6pOH2WzmBhG5E8hT1dfquN0zwBRV/a5xIjPG/1h9ajh2BWUAXuLX7eTuWmuVyZjfsPrUQOwKyhhjjE+yKyhjjDE+yRKUMcYYn2QJyhhjjE+yBGWMMcYnWYIyxhjjkyxBGWOM8UmWoIwxxvgkS1DGGGN8kiUoY4wxPskvlttISEjQ5ORkb4dhDADLli07oKqtvB1HfVl9Mr6mujrVaAlKRN4Azgb2q2pf13stgY+AZJwp5S9V1eza9pWcnMzSpUsbK1Rj6kREdnroezoC7+AswVAJTFLVZ44ocyLwJbDd9dbnqvpoTfu1+mR8TXV1qjGb+N4Cxh7x3p+B71S1O/Cd67UxpmrlwF2q2gsYCdwiIr2rKDdfVQe6HjUmJ2P8SaNdQanqPBFJPuLt84ATXc/fBuYA99X3O9bM+x+Fq/6HBoahQWEQFAbB4UhwOAEh4UhIJIGhkQSFRhIUHkVwWBQhETGEREQTHhVLeHg4IYEBHLaEszE+w7UEeZrreZ6IbADaA+u9GpgxblJVvlyZyub0PO4d27PO23v6HlQbV6VDVdNEpHV1BUVkIjARICkpqcoyhfu20C1zDqFaSiilBEtFnYIp1UCyCadQIiiScIoCIikOiKI0KIrS4GjKg2OpCI2lMiwOiWhBQGRLgqPiCY1pRXhMK2IjQ4kNDyE2PJiQIOtvYhqP62RvELCoio9HicgqIBW4W1XXVbF9rfXJmIa0dm8OD09Zx9Kd2QzoGEdJeQWhQYF12kejLrfhqlRfH3YP6qCqxh32ebaqtqhtP0OHDlV32swryssoLsynpKiA0uJ8SosKKCsuoLwon/KSAsqL86gsyaeypAAtyUdKnUdgWT6BZXmElBcQUp5PeGU+EZX5RGoBAVT9+6lUIZsoMjWGTI0lOyCWvKAECkPiKQlvTUVkGzQ6keAWHYiJbUF8ZAgJ0aG0jg4lISqUsOC6/UMZ3yEiy1R1qAe/LwqYC/xDVT8/4rMYoFJV80VkHPCMqwm9Wu7WJ2Pq40B+CU/O2MRHS3cTHxnCvWf05OIhHQgIqL6lqro65ekrqHQRSXRdPSUC+xty54FBwUTGtCAyptac557KSijJRYuyKcnLojgng+K8TMryMqjIz0ALDhBakEGn4gP0KNlNZNlKwoqKoAjI+mU3BzWSVE1gryawShPYra3JDEmkMLITlS2SSYiNITEujHax4bSLC6dDi3AS48LqfLZhmh4RCQY+A94/MjkBqGruYc+niciLIpKgqgc8GacxZRWVvLNgJ0/P2kxRaQU3jO7Mbad2JyYsuN779HSCmgJcAzzh+vmlh7+/bgICIDwOCY8jrGVnwtzZpiQf8vZBXirkplF+cA/BmbvokL2bpNzdhOavJ7iiCBTIh8p8Yd/uBLZUJLJV2zFd27GpsgNb6EB4dAJJ8RF0ahlBckIknRMi6dIqkuT4SLsCawbEuTn6OrBBVf9bTZm2QLqqqogMx+n4lOnBMI1h3uYMHv16PSn78xnToxUPnt2bbq2jjnq/jdnNfDJOh4gEEdkDPISTmD4WkRuAXcAljfX9XhMaBaHdIKEb4PyCf/VLVoXCLMjeAVlbCcjaRrsDW2h7YAvHZ84joKzw56LZla3YmtmZlelJLC7uyDuVXUmnJSKQ1DKC7q2j6NY6ml6J0fRsG0OXVpEEB9q9sCZkNHAVsEZEVrreewBIAlDVl4GLgT+ISDnOtft4tWWyjYfsyizk71PXM3N9Op3iI3j16qGc2qt1g3U884sl35tNm7kq5OyG/Rth/3pIXwv71sCBzaCVABSFtWZvRG/WBPRkbnFXvs1uS2GFczUVEhhAz8Ro+rSLpV/7WAZ2jKNHmyiCLGk1KE/fg2pozaY+mUZTUFLOi3NSeHX+doIChD+e3I0bjutc79sSvnIPytREBOKSnEeP0395v7TQSVZ7lxOeupxuuxfT7cAcLgA0PJzCNkPYGTOExdKPmQdjmLo6lcmLdwEQERJI/w6xDEtuybDklgzu1IKoUPtnN8bUnaoyZVUqj0/byL7cYi4Y1J77xvakbaxbN0DqzP6n8gchEdBxuPM4JC8ddi9Cdv5E5I759N7wDL2Ba8Pi0D4nk5l4PEuCh7EoPYBlO7N5YXYKlQqBAUL/DrGM7prAsd3iGdqppXWRN8bUau3eHB75ah1LdmTTr30sL0wYxJBOLRv1O62Jr6koOADb50LKd5AyC/LTQQIgaRT0PIuCbmezLDuCRdsz+WlrJqv35FBRqUSGBDK6WwIn9WzNKT1b0zqmcc6EmhJr4jPNSWZ+CU9+u5kPl+yiZUQI9449hkuGdKyx23hdVVenLEE1RZWVsG81bJoGG6c6zYPgJKu+F0Hfi8gNiGbh1kzmbs5g9sb9pOYUAzAoKY4z+rRlXN9EkuIjvHgQvssSlGkOyioqeW/hTp6auZnC0gquHpXM7ad2Jza8/t3Gq2MJqjnL3ArrPoc1n0HGBggIhp7jYOCV0O0UVALYlJ7HzHXpzFi/j7V7naE1AzvGcc6AdpwzIJHW0XZldYglKNPU/bDlAI98tY4t+/M5vnsCD53Tm26toxvt+yxBGce+NbDyA1j9ERRmQmwSDLseBl0NkfEA7M4qZOqaNL5alcq61FwCA4QTe7TikqEdObln62Z/z8oSlGmqdmcV8tjU9cxYl05Sywj+elYvTuvdptHnK7UEZX6tvNRpAlzyGuyYD4GhMGA8HHvbz2O4AFL25/Ppsj18vnwP+/NKaBUdyuXDk7hieFKj9dzxdZagTFNTWFrOS3O28sq8bQTKL93GPTUhgCUoU7309bD4FVg5GSpKodc5MOZuSBzwc5Hyikrmbs7gvYU7mbM5gwARzuzblpvGdKVfh1gvBu95lqBMU6GqfLU6jcenbSAtp5jzB7bjvjN7khgb7tE4LEGZ2uXvh0UvO1dVxTnQ61w46QFo3etXxXZmFvDewp18uHg3eSXlHNs1nptP7MbobvHNYukSS1CmKViXmsMjU9azeEcWfdvH8PA5fRia3LjdxqtjCcq4rzgHFrwIC16A0nyn6e+UByGm3a+K5RaXMXnRLt74cTvpuSUMS27Bnaf14NiuCV4K3DMsQRl/llVQypPfbuLDxbuIiwjh3jOO4ZKhHQlswG7jdWUJytRdYRb88JRzVRUQBKNvh2NvhZDIXxUrKa/goyW7eWF2Cum5JRzbNZ4HxvWib/um2fRnCcr4o3JXt/H/ztxMQWkFV4/qxB2n9miUbuN1ZQnK1F/2Dpj5EKz/wun1d9Z/fj0Vk0txWQUfLNrFc99v4WBRGRcMas89Zxzj8fbsxmYJyvibH1OcbuOb051u4w+e3ZvubRqv23hdVVenmnd/YeOeFslw6dtw7TQIDocPLoGPr3GWFTlMWHAg1x/Xmbn3nsTEMV34enUaJz85l1fmbqWsotI7sRvTjO3OKuT37y5jwmuLKCqrYNJVQ3jn+uE+lZxqYldQpm7KS+GnZ2Du/znJ6uz/OrNTVGF3ViGPfu1Mxd+jTRSPnd+P4Z29cxO2IdkVlPF1RaUVvDQnhVfmbSPAC93G68quoEzDCAqBMffAH36C+K7w6fXw6Q3O/aojdGzprA/z2tVDKSyt4NJXFvDwlHUUlpZ7IXBjmj5V5atVqZzynzk8+30KZ/Rpy/d3n8AtJ3Xz2eRUE5vN3NRPQje4/lunE8XcJ2DXQrjkLeg47DdFT+3dhmO7xfPvbzbx1k87mL1pP09eMoBhXurSakxTtD41l4e/Wsfi7Vn0aRfDM5cP8vs6ZldQpv4Cg+CEe+CGmRAQAG+OdbqmV9FsHBESxMPn9uHDiSNRhcteWcBTMzdTbvemjDkq2QWl/PWLNZz93HxS9ufzzwv6MeWPx/l9cgJLUKYhtB8MN82D7mfAjAfg46uhtKDKoiO7xDP99uO5cHAHnvluC1e8uojUg0UeDtgY/1deUcnbP+3gxCfnMHnxbq45NpnZd53IFSOSvDqmqSFZgjINI7wFjH8fTvs7bPwa3jgDcvZUWTQyNIgnLxnAU5cNYG1qDuOenc+PKQc8HLAx/uunlAOc9ewPPDRlHf3axzL99uN56Jw+xEZ4f0xTQ7IEZRqOCIy+DS7/CLJ2wKSTYE/1vcUuGNSBr289jtbRoVz1+iJem78Nf+hVaoy37M4q5A/vLeOK1xZRUFrOK1cN4d0bhtPDT7qN15UlKNPwepwON850uqG/dTZsnlFt0S6tovj85tGc3rstj03dwJ0fraS4rMKDwfouEekoIrNFZIOIrBOR26soIyLyrIikiMhqERnsjVhN4yoqreC/Mzdz6n/nMmdTBnef3oNZfzqBM/q0bdLzX1qCMo2jdS+48Tto1QMmXw6rPqy2aFRoEC9OGMxdp/Xgi5WpXP3GYg4WlnowWJ9VDtylqr2AkcAtItL7iDJnAt1dj4nAS54N0TQmVeXr1a5u499t4XRXt/E/ntzdL7uN15UlKNN4olrBNV9D8mj4303OBLTVCAgQbj2lO89ePoiVuw5y0Us/sTur0IPB+h5VTVPV5a7necAGoP0Rxc4D3lHHQiBORBI9HKppBOtTcxk/aSF//GAFcREhfHzTKJ67fFCTmzqsJpagTOMKi4EJn0Lv82DG/fDjMzUWP3dAO969YTgZeSVc8OJPbEjL9VCgvk1EkoFBwKIjPmoP7D7s9R5+m8QQkYkislRElmZkZDRWmKYBZBeU8rcv1nL2c/PZnJ7HPy7oy1e3HtckZmGpq1oH6orIUOB4oB1QBKwFZqnqb6cOMKYqQaFw0Rsgv4OZDzrvjf7N7ZSfjegSz+c3H8uVry1m/KSFvHP9cAZ0jPNQsL5HRKKAz4A7VPXIjF3VDYjf9DRR1UnAJHCmOmrwIM1RK6+o5IPFu/jPt5vJLynn6lHJ3HlqjybXM68uqr2CEpFrRWQ5cD8QDmwC9gPHATNF5G0RSfJMmMbvBQbBha9CnwudJPXjszUW79Y6mk9+P4qY8CAmvLaIJTua5/mQiATjJKf3VfXzKorsAToe9roDkOqJ2EzDWbA1k7Of+4EHv1xH3/YxTLvteB4+t+l1G6+rmq6gIoHRqlrlKEoRGYhzY3ZXYwRmmqBDSQpg5t8gLBaGXFNt8Y4tI/j4plFMeG0RV7++mLeuG8aILvEeCtb7xOme9TqwQVX/W02xKcAfReRDYASQo6ppnorRHJ092YX8c9oGpq3ZR4cW4bx85eAm3zOvLqpNUKr6QnWficgdqvp044RkmrTAILhwEpTkwdd3QERL6HVOtcUTY8P5aOIoxk9awPVvLeG9G0cwKKmFBwP2qtHAVcAaEVnpeu8BIAlAVV8GpgHjgBSgELjOC3GaOioqreDluVt5ee5WROBPp/Vg4pguzaJnXl3Ua7kNEdmlqh5r3rPlAZqg0gJ45zxIWw1Xfgadj6+xeHpuMZe+soDsglI++N1Ir67WW5/lNkSkNU7COfxe7lJV9fhkhFafvEdVmbZmH/+ctoG9B4s4u38iD4zrRbu45tMzryoNvdyGXX+aoxMSCVd8DC07w4dXwP6NNRZvExPG+zeOIDosmKteX0TK/nwPBXp0ROQkEZkBTMUZs5QI9Ab+inNl9IiIxHgzRuMZG9JyufzVhdzywXJiwoP5aOJInr9icLNPTjWpb4KyXkDm6EW0dLqgB4XB5MugILPG4h1aRPD+jSMIDBCufXMx+/OKPRToURkH/E5Vh6nqRFX9q6rerarnAgOAFcBp3g3RNKaDhaU8+OVaznp2Ppv25fHY+X35+tbjmtX91PqqqRdfnojkHvbz0CMPp5nCmKMX1xHGfwC5afDxVc6KvTVITojkjWuHkVVQynVvLiG/xLcXP1TVe1S1uo5E56nqF6r6mUeDMh5RXlHJuwuc2cbfX7SLq0clM/vuE7lyZKcmM9t4Y6s2QalqtKrGHPYz5rDXttChaTgdh8H5L8LOH2HqnVWuJ3W4/h3ieGHCYDbuy+MP7y2jzH/XlHrK2wGYxrFwm9Nt/G9frqN34i/dxuMiQrwdml+p6QpqqYg8IyJjRSTMk0GZZqjfxTDmXljxHix7s9biJx3Tmscv6Mf8LQf4+9frPRBgo7DT6CZm78EibvlgOeMnLSSvuJyXJgzm/RtHcEzbpjnbeGOr6UpoJM6g3LHAIyKSCcwApqvqZk8EZ5qZE++H1OUw/T5IHOgshFiDS4d1ZGtGPq/M20bPtjFcMcLvxo3bvdwmorjsl27jYN3GG0pN46DKgTmuB64JKM8EHhOR7sACVb3ZAzGa5iIgwBnI+8oY+PgauGmu05GiBveO7cmm9Dwe/HItXVtF+tyNZxFZQ9WJSIA2Hg7HNDBVZfraffxjqtNt/CxXt/H21jOvQbh9L8k1Ov0N4A0RCQBG1fdLRWQHkAdUAOV1HVNimrCIlnDJ286KvJ9PdLqiB1Tf2TQwQHhm/CAuePFH/vD+cr669Thf+8/hbG8HYBrHxn25PDJlPQu2ZdKzbTQfThzJSB87QfJ37kwW+xW/PQPMAZa6BlfVt6/vSapq63yb3+owBMY+DtPuhkUvw6iaL9Rjw4N57eqhnPv8j9zy/nI+vmkUIUE+M1H/Lq1lNLyISG1ljO84WFjKUzM38+7CncSEB/P38/ty+bCOBAX6zN9ck+HOb3QbkA+86nrkAulAD9drYxresBvhmHEw6yHYt6bW4l1aRfHvi/uzcvdBnphe86BfD5stIrceObGyiISIyMki8jZQ/YSExmdUVCrvLtzJSU/O4d2FO7lqZCfm3H0iV43sZMmpkbjTxDdIVccc9vorEZmnqmNEZF09v1eBb0VEgVdcSwH8iohMxFkhlKQkv7v5bY6WCJz7HLx0LHx2I0yc4ywhX4Nx/RK59thk3vhxO8M7t2BsX59Yt28scD0wWUQ6AweBMCAQ+BZ4SlVX1rC98QGLtmXy8Ffr2ZCWy6gu8Tx0bm96trUJQBqbOwmqlYgkHRps6DoTTHB9Vt91uUeraqprfrKZIrJRVecdXsDWrzFEJjjjo967yFmiY9z/1brJA+N6sWL3Qe75ZDW9EmPoFB/pgUCr52oCfxF40bV0RgJQpKoHvRqYcUvqwSL+OW0DX69Oo31cOC9OGMyZfW22cU9x57r0LuAHEZktInOA+cA9IhIJvF2fL1XVVNfP/cD/gOH12Y9pBrqdCiNvhsWTYOv3tRYPCQrghSsGIQJ3frSSch8axKuqZa5l3C05+bjisgqembWFk/8zh5nr07nj1O7M+tMJjOuXaMnJg2pNUKo6DWfdpztcj2NUdaqqFtRnyQ0RiRSR6EPPgdNxZnY2pmqnPAjx3WHKbc4yHbXo0CKCxy7ox/JdB3l+dooHAjRNhaoyfU0ap/xnLk/N2swpPdvw3V0ncMepPQgPsTFNnlbTTBLHHXquqiWqukpVVx7qtSciMSLStx7f2QbnimwVsBiYqqrf1GM/prkIDofzXoCcPTDzIbc2OXdAOy4Y1J7nvk9h+a7sRg7QNAWb9uUx4bVF/OH95USHBTH5dyN5YcJgOrSI8HZozVZN96AuEpF/A98Ay4AMnJu73YCTgE44zX91oqrbcGZxNsZ9SSOcpr6FL0Dv86DLCbVu8sh5fVi8PYs7P1rJ1NuOJyrU81NIisjzwAeq+pPHv9y45VC38fcW7SIqNIhHz+vDFcOTrGeeD6hpstg7gbOANOAS4O/An3Ca+15R1TGqusQjURoDcPJfoWUXmHKrs+BhLWLCgnl6/EB2ZxXy+LQNHgiwSluA/4jIDhH5l8j/t3fn8VVV1wLHf4tMjIJhDEMAmSRhLIEEBEQQRBRqW62AT3xWS8Uq2qp9rdoCFp9ofW3pQ6jW4lQcanFAREGoDFUikygQ5jEMJiBzCUiS1T/2iQkQkktI7rn3Zn0/n/tJ7j3n3OwN2Vl3n7PO2tLFr4aYM+XlKzM+K0wbH9kjkYUP9mNUzxYWnMrT0b2QuaxMh5b4kVJVD1F4/5Mx/oqt7k71vXAtLHwCBk0s9ZDuLeK5o3dL/rJkO9d1TKBX63qlHlOeVHUyMFlEmgPDgRe84suvAa9bXUt/LNt+kPGz1pGx7yipLeMZPyyZ9gmWNn7R8vNh/wbITIdd6bBrKRzeBbUaw88z3O0jF8CWzTDhpXkv6HorLJ0KnW6GRh1LPeSBQe1YsD6bECwa2gAAEr1JREFUX8z8krn396WGD6f6VHUn8CTwpIh0xZUNG4e7H8oEiaWNl7PTJ12B511eQMpMh5NH3LYaDaB5T0gdA4lpZXp7C1Am/Ax8DDbOgdk/gx/NK7FWH0DVmCieurETNz27lKc+3MCE75Ylt+fiePdADcbNogYAi4AJQW9IJXXydB5/WbyNZxZuQRXuG9CGu65sZZl5F+rEQcj8zM2MdqXD3s8hz7sdtl5baD/MfYhMTINLW17wjOlsFqBM+KkeD4Meh3fuglUvQsqPSj0kpUU8t/dqyfRPtjOkY0LQqp6LyEBgBO567jLgdWC0qpZ6EU1EpuOKzWar6jlRVUT6Ae8C272X3lLVx8qp6RFBVZm77ismvr+e3YdyGNKxEQ8PaW+ZeYFQhUPbYVeRgHRgo9tWJcYth5N6lwtGzdKgRvmPqfMGKBHpDmSq6lfe81HAD4CdwHhVPVjurTEmUJ2Hw+oZMH88XH491GxQ6iEPXdOO+euz+NXba5gztk+w1up5GHgVeLAMY+ZFYArwcgn7LFFVq5hejE1Zx5jw3jo+2fI17RrW4tUfp9KrVXCvQYaVvFzIWlN47WhXOhzPctviakOzHtD5ZkjsCY27llp6rDyUNIN6FrgaQET6ApOAe4EuuBJEN1Z464w5HxG4/g8wtacLUjdMLfWQarFRTLyhA6OmL2Pqwq38fGDbCm+mql51EccuFpEW5deayuHIidP8Yb6rNm5p4yU4dRx2Ly8MSLtXwGlvYl8nEVpe6WZHiT2h/uWlnkqvCCUFqKgin/huBp5T1ZnATBGx4pbGf/XauKU4PpnsTvM1LX1Zsb5t63NDl8ZMW7iFoZ0SaNMw7Jfi7und9L4XN0srtoBzZSi+nJevvLE8k6fnbeTwiW8YmZrIAwPbcWmNWL+bFhqO7iuSXZfuVgnQPJAq0DAZut4CzVJdQKrdxO/WAqUEKBGJ9lbWHYD3yx3AccYET9+H4Is3YM5DcOeCgD7lPXp9Egs37efht9fwxuieVKkSthlcq4DmqnpcRIYA7+DuUzxHpBdfXr7DpY2v23uUHi3jGT80maTGlThtPD/fXS/atbTwGtLhnW5bTHVo0g36POBmSE27Q9XQ/LcqKdC8BiwSkQNADq5ILCLSGrdgoTH+i6vlsvreHg2r/wbfGVXqIfVqxvHwkPb84h9f8vryTEamhueMQlWPFvl+johMFZF6lWkh0H1HcnhizgZmfbGXxrWrMmVkV66rjAVdT5+EfasLrx3tSoeTXk3iGg1cJZYeo13ad6NOEBXjb3sDdN4ApaqPi8gCIAGYV2TFzyq4a1HGhIZOP4QV02H+BJfmWq1OqYfc1K0pM1fu5skPN3BNckPq1owLQkPLl4g0ArJUVUWkB25sfu1zs4Li5Ok8nl+yjWc+3kq+KmMHtGFMZUobP3HQVWf4Nt17VWG6d93W0H6oO1WXmOaqr4RpwC6tkkR6Ma/Zne8mtIjAkKfg2Sth0VMw+H8DOESYeEMHrp28hCc+2MDTN4VeeUgReQ3oB9QTkd24G3tjAFT1z7hEpTEikos7yzE80peOd2njWTw+J4PMgzlc28GljTeLj+C0cVU4tOPM+4/2e6tGV4lxGXWpP3EBqVmqW0ctQti1JBMZEjrDd26FZc+6hIl6rUs9pE3DWtzZ5zL+vGgrP0xpRo+W8UFoaOBUdUQp26fg0tArhc1Zx5jwXgb/2nLApY3fmRr00lVBkZcLWWvPSvf+ym0rSPfueJObHTXpFpR0b79YgDKRo/+vYe3bMO9RGPl6QIeMHdCa977Yy6/fWcvssb2JsVTkkHMk5zR/nL+Jl5fupEZsFBOGJXNLagSljZ86DntWnJnu/c1xt612IrTsUyTdu70v6d5+sQBlIkfNBtD3QZg/zq2+26p/qYdUj41m3NAkRr+ykhc+2c7ovq2C0FATiLx85e8rMvndXJc2PqJHIg8Makd8uKeNH/vqzNp1+7506d4INOwAnUd4ASkNajf1u7W+sgBlIkvaGFj5Asx9BH6yBKJK/xUflNyIAZc34I/zNzO0c2MSakfuKZNwsWLHQcYVpI23iGfcsCSSG9f2u1kXLj8fDmw6s7r3oR1uW3Q1d+9en58XSfcOwz5WIAtQJrJEx8HA38Lfb3V1+rrfGdBh44clc/XvF/Hb2RlMvaVbxbbRnNe+IzlM+mAD767eS0Ltqvz/iK5c3ymM0sZzT8HeIunememQ463oXL2eC0Td74TEXpAQPunefrEAZSJP+6HQvDf883HocGNAaefN4qtzz1Wt+b+PNrF40376tq0fhIaaAidP5/HXf21nyj+3kKfK2P6tuatfK6rHhvifqJxDZ6Z771kFeafctvhWcPl1XnZdGtRtFbbp3n4J8f99Y8pABAY/Ac/2DTjtHGD0lZfx1ud7GDdrHR/e34e46EpyT42PVJV5GVlMfN+ljQ9ObsQj14Vo2riqq8ZQUJkh8zPIznDbqkRDQhfo8ePCdO+a9iHnYlmAMpEpodMFp53HRbsMsVHTl/Hcom3cO6DYqkGmnGzOOsZjszNYsvkAbRvWZMadqVwRSmnj+Xnnpnsf2+e2xV3iglCH77vZUZNubsVnU64sQJnI9W3a+SMw8o2ADunbtj6jejanZf0aFdy4yutIzmkmz9/MS0t3UCM2ivFDk/ivtOb+p41/82+X4n1Guvcxt+2SptD8isLsugZJUMVm2BXNApSJXDUbQN8H3HIcW+ZD66sDOuwxH1bcrQzy8pU3vbTxgwVp4wPb+ldm6liWl13nnbLb90WRdO9kV0KreS83U6rTzJ82VnIWoExkS7sbVr0MH/wPjPnUZfmZoFux4yDj31vH2j1H6d7iUl4a2oMOTYKYUq3q0r13FU339hYijq4KTVKg988K070DSKwxFc8ClIls0XFw7e9gxg9g6RS3xIAJmq+OnGTSB+t5Z/VeGl1SlT+N6MrQYKSN555yM6Ki1b1zvOXtqtd114263+G+JnSG6DC/+TdCWYAyka/N1W5Z+EW/czXM6oTn8hrhpCBt/JmPt5Cbr9zbvzVjKjJtPOfwWeneK89M9243xC05kdjL0r3DiAUoUzkMngRTusOHv4LhM/xuTcRSVeavz+a3szPYdfAE1yQ35NHrkso3bVwVjmQWya4rSPdWl+7dqJNL926W6k7Z1WxQfj/bBJUFKFM51GkGVz4ECx6DTXOh7TV+tyjibMl21caXbD5AmwY1+dsdqfRuUw5p4/l5kLXurHTvvW5bbC1X3Tv5Bq+6d4qle0cQC1Cm8uh5L6z5B8y6F8YshRp1/W5RRDh60ksb/3QH1WKjGOeljZe5Mvw3/3an6AoCUubyIuneTdyqsAU3wzZMtnTvCGYBylQe0bHw/efgL/1h9n3ww1fsWsRFyM9X3lyZyVMfurTx4d0TeXBQGdLGj2d7deuKpHvn5wLi7jfqdJO7dpSYZunelYwFKFO5NOoI/R+Fj34Dq1+Frrf43aKwtHLnQcbPymDNniOkNL+Ul4YFmDauCgc2F6nunQ4Ht7ptUXGuuvcV97nsumY9LN27krMAZSqfnvfApnnu3qjmvSC+pd8tChtZR08y6YMNvP35HhpdUpXJw7swrHPj86eN535zZrp3Zjqc+NptqxbvZkXdbnOn7BI6231q5gwWoEzlUyUKvjcNpvWGGTfC7R9YplcpTuUWVhvPzVd+elUr7u7XmhpxZ/0JyTkMu5efme6de9Jti78M2g72sut6Qr02dorVlMgClKmc6iS6+nyvfM89bnsPqsf73aqQo6osWJ/NxPcz2PH1CQYlubTxxLpeptzhIunemZ+5bDsUJMoV7E25w91/1CwNajX0tS8m/FiAMpVX854w4lV49WY3kxr1LsTV8rtV3xKR6cD1QLaqnlMgUNx5tcnAEOAE8N+quqq8fv6W7OM8NjuDxZv207pBTV65vRt9au+HrTPgY2+GdHSP2zm2pisR1O+XbnbUNAVireCuuTgWoEzl1qo/3PQivHErTLsCrvu9qzwRGl4EpgAvn2f7tUAb75EKTPO+XpSjJ0/zp/mbef3TjXSP3c7MpAN0lQ1UeWs5nDrqdqrVuLCyd2IaNEiGKPtzYsqXL79RIjIY98kvCnheVSf50Q5jALfq6W2z4L37Xc2+5O+7BQ9rNfK1Waq6WERalLDLd4GXVVWBdBGpIyIJqrqvLD8vP+coS+fPZPuqBVyfl8GvYncQRR5sA+q3h443Ft5/VCfRrh+ZChf0ACUiUcAzwEBgN7BcRGapakaw22LMt1r0hjGfwCeTYfHTEBXj7pkKbU2AzCLPd3uvnROgRGQ0MBogMbH4WoQbN63nipX3050YTjfuSlTrYV5A6gHVLi3/1htTCj9mUD2ALaq6DUBEXsd9ErQAZfwVHQdX/gI6/ABiqvndmkAUN4XR4nZU1eeA5wBSUlKK3ad9x+6sy3mTpG59iY2pWn6tNKaM/AhQxX3qO+e8eSCf+IypEHVb+d2CQO0GipZWaArsLfO7ValCctqgi22TMeXGjzWWA/rUp6rPqWqKqqbUr18/CM0yJuzMAkaJkwYcKev1J2NCkR8zqPL91GdMhBKR14B+QD0R2Q2MA2IAVPXPwBxcivkWXJr57f601JiKIS4BKIg/UCQa2AQMAPYAy4GRqrquhGP2AzvPs7kecKC82+kD60doKakfzVU1bKf1pYwniIz/w0joA1SefhQ7poI+g1LVXBG5B5iLSzOfXlJw8o457x8DEVmhqinl3Mygs36ElkjpR3FKC66R0PdI6ANYP3y5D0pV5+BOTxhjjDHF8iNJwhhjjClVJASokL+bMkDWj9ASKf0oi0joeyT0ASp5P4KeJGGMMcYEIhJmUMYYYyKQBShjjDEhKawDlIgMFpGNIrJFRH7pd3sCJSLTRSRbRNYWeS1eRD4Skc3e15CuzikizUTkYxFZLyLrROQ+7/Vw60dVEVkmIl94/Zjgvd5SRD7z+vGGiMT63daKZuPJXzamzhW2AapIVfRrgSRghIgk+duqgL0IDD7rtV8CC1S1DbDAex7KcoEHVLU9kAb81Pv3D7d+nAL6q2pnoAsw2Csb9CTwB68fh4A7fGxjhbPxFBJsTJ0lbAMURaqiq+o3QEFV9JCnqouBg2e9/F3gJe/7l4AbgtqoC6Sq+wpWb1XVY8B6XCHgcOuHqupx72mM91CgP/AP7/WQ70c5sPHkMxtT5wrnAHW+tXDCVcOCQp/e1wY+tydg3qJ6XYHPCMN+iEiUiKwGsoGPgK3AYVXN9XYJ99+tQNh4CiE2ppxwDlABr4VjKo6I1ARmAver6lG/21MWqpqnql1whYt7AO2L2y24rQo6G08hwsZUoXAOUJFWFT1LRBIAvK/ZPrenVCISgxtIM1T1Le/lsOtHAVU9DCzEnf+v4xU2hvD/3QqEjacQYGPqTOEcoJYDbbzMkFhgOG59nHA1C7jN+/424F0f21IqERHgr8B6Vf19kU3h1o/6IlLH+74acDXu3P/HwI3ebiHfj3Jg48lnNqaKoaph+8CthbMJd37zEb/bcwHtfg3YB5zGfXK9A6iLy9DZ7H2N97udpfShN26K/iWw2nsMCcN+dAI+9/qxFviN9/plwDLcWktvAnF+tzUI/xY2nvzth42psx5W6sgYY0xICudTfMYYYyKYBShjjDEhyQKUMcaYkGQByhhjTEiyAGWMMSYkWYAKMSJSV0RWe4+vRGRPkeefVtDP7Coiz3vfjxeRBwM8bn6oV1Y2lZuNp/AWXfouJphU9WtcBWBEZDxwXFWfruAf+zAwMdCdvRsKBXgFuBt4vILaZcxFsfEU3mwGFUZE5Lj3tZ+ILBKRv4vIJhGZJCK3eGuwrBGRVt5+9UVkpogs9x5XFPOetYBOqvpFkZeTRGShiGwTkbHefi28dWqmAqtwZXFmASMqut/GVAQbT6HPZlDhqzOuAONBYBvwvKr2ELfI2b3A/cBk3Por/xKRRGAu5xZtTMHd7V3U5cBVQC1go4hM815vB9yuqncX7CgicSJS1/ukaky4svEUgixAha/l6pXgF5GtwDzv9TW4wQCuBlaSO4MAwCUiUkvdWjMFEoD9Z733+6p6CjglItlAQ+/1naqafta+2UBjwAaUCWc2nkKQBajwdarI9/lFnudT+P9aBeipqjklvE8OULWE984r8n7/Lub4qt57GBPObDyFILsGFdnmAfcUPBGRLsXssx5oXZY39y7uNgJ2lOV4Y8KMjacgswAV2cYCKSLypYhkAHedvYOqbgBqexd3L1Q3IF0LV8k0JpLZeAoyq2ZuEJGfAcdU9fkLPG4yMEtVF1RMy4wJPzaeyo/NoAzANM48Tx6otTaYjDmHjadyYjMoY4wxIclmUMYYY0KSBShjjDEhyQKUMcaYkGQByhhjTEiyAGWMMSYk/Qf5waP3lkrOWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = {}\n",
    "\n",
    "fig, [[ax['X'], ax['P']], [ax['S'], ax['V']]] = plt.subplots(2, 2)\n",
    "\n",
    "for F, out in zip(Fs, results):\n",
    "    var = {name: y for name, y in zip(names, out.y)}\n",
    "    for name in names:\n",
    "        ax[name].plot(out.t, var[name])\n",
    "        ax[name].set_ylabel(f'{name} ({units[name]})')\n",
    "        ax[name].set_xlabel('Time (hr)')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}